library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)
library(ggplot2)
library(rstatix)
## 
## Attaching package: 'rstatix'
## The following object is masked from 'package:stats':
## 
##     filter
library(ggpubr)
library(dbscan)
## 
## Attaching package: 'dbscan'
## The following object is masked from 'package:stats':
## 
##     as.dendrogram
library(corrplot)
## corrplot 0.94 loaded
library(corrr)
library(factoextra)
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
theme_set(theme_minimal())

1

Сделайте копию датасета, в которой удалите колонки с количеством пропусков больше 100, а затем удалите все строки с пропусками.

df_full <- readRDS("~/BioStat/Прод_виз_данных/very_low_birthweight.RDS")

df_cleaned1 <- df_full %>%
  select(where(~ sum(is.na(.)) <= 100)) %>%
  filter(if_all(everything(), ~ !is.na(.)))
summary(df_full)
##      birth            exit          hospstay            lowph      
##  Min.   :81.51   Min.   :68.53   Min.   :-6574.00   Min.   :6.530  
##  1st Qu.:83.52   1st Qu.:83.58   1st Qu.:   16.00   1st Qu.:7.130  
##  Median :84.90   Median :84.96   Median :   37.00   Median :7.210  
##  Mean   :84.75   Mean   :84.84   Mean   :   40.36   Mean   :7.202  
##  3rd Qu.:86.07   3rd Qu.:86.17   3rd Qu.:   62.00   3rd Qu.:7.310  
##  Max.   :87.48   Max.   :96.87   Max.   : 3668.00   Max.   :7.550  
##  NA's   :21      NA's   :31      NA's   :31         NA's   :62     
##      pltct                    race          bwt            gest      
##  Min.   : 16.0   black          :369   Min.   : 400   Min.   :22.00  
##  1st Qu.:143.0   native American: 16   1st Qu.: 900   1st Qu.:27.00  
##  Median :202.0   oriental       :  4   Median :1120   Median :29.00  
##  Mean   :201.6   white          :257   Mean   :1094   Mean   :28.87  
##  3rd Qu.:252.0   NA's           : 25   3rd Qu.:1310   3rd Qu.:31.00  
##  Max.   :571.0                         Max.   :1580   Max.   :40.00  
##  NA's   :70                            NA's   :2      NA's   :4      
##           inout          twn              lol             magsulf      
##  born at Duke:547   Min.   :0.0000   Min.   :  0.000   Min.   :0.0000  
##  transported :121   1st Qu.:0.0000   1st Qu.:  0.000   1st Qu.:0.0000  
##  NA's        :  3   Median :0.0000   Median :  3.500   Median :0.0000  
##                     Mean   :0.2074   Mean   :  8.438   Mean   :0.1344  
##                     3rd Qu.:0.0000   3rd Qu.:  9.000   3rd Qu.:0.0000  
##                     Max.   :1.0000   Max.   :192.000   Max.   :1.0000  
##                     NA's   :20       NA's   :381       NA's   :247     
##       meth             toc              delivery        apg1      
##  Min.   :0.0000   Min.   :0.0000   abdominal:314   Min.   :0.000  
##  1st Qu.:0.0000   1st Qu.:0.0000   vaginal  :335   1st Qu.:2.000  
##  Median :0.0000   Median :0.0000   NA's     : 22   Median :5.000  
##  Mean   :0.4372   Mean   :0.2248                   Mean   :4.903  
##  3rd Qu.:1.0000   3rd Qu.:0.0000                   3rd Qu.:7.000  
##  Max.   :1.0000   Max.   :1.0000                   Max.   :9.000  
##  NA's   :106      NA's   :106                      NA's   :34     
##       vent            pneumo            pda              cld        
##  Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000  
##  1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000  
##  Median :1.0000   Median :0.0000   Median :0.0000   Median :0.0000  
##  Mean   :0.5803   Mean   :0.1969   Mean   :0.2087   Mean   :0.2694  
##  3rd Qu.:1.0000   3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.:1.0000  
##  Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1.0000  
##  NA's   :30       NA's   :26       NA's   :29       NA's   :66      
##        pvh            ivh            ipe           year           sex     
##  absent  :360   absent  :442   absent  :472   Min.   :81.51   female:320  
##  definite:125   definite: 75   definite: 38   1st Qu.:83.52   male  :330  
##  possible: 41   possible: 10   possible: 17   Median :84.91   NA's  : 21  
##  NA's    :145   NA's    :144   NA's    :144   Mean   :84.76               
##                                               3rd Qu.:86.07               
##                                               Max.   :87.48               
##                                               NA's   :21                  
##       dead       
##  Min.   :0.0000  
##  1st Qu.:0.0000  
##  Median :0.0000  
##  Mean   :0.2146  
##  3rd Qu.:0.0000  
##  Max.   :1.0000  
## 
str(df_full)
## 'data.frame':    671 obs. of  26 variables:
##  $ birth   : num  81.5 81.5 81.6 81.6 81.6 ...
##  $ exit    : num  81.6 81.5 81.6 81.7 81.6 ...
##  $ hospstay: int  34 9 -2 40 2 62 32 NA NA 28 ...
##  $ lowph   : num  NA 7.25 7.06 7.25 6.97 ...
##  $ pltct   : int  100 244 114 182 54 NA 282 NA NA 153 ...
##  $ race    : Factor w/ 4 levels "black","native American",..: 4 4 1 1 1 4 1 NA NA 1 ...
##  $ bwt     : int  1250 1370 620 1480 925 940 1255 600 700 1350 ...
##  $ gest    : num  35 32 23 32 28 28 29.5 26 24 34 ...
##  $ inout   : Factor w/ 2 levels "born at Duke",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ twn     : int  0 0 0 0 0 0 0 NA NA 0 ...
##  $ lol     : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ magsulf : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ meth    : int  0 1 0 1 0 1 1 NA NA 1 ...
##  $ toc     : int  0 0 1 0 0 0 0 NA NA 0 ...
##  $ delivery: Factor w/ 2 levels "abdominal","vaginal": 1 1 2 2 1 1 2 NA NA 1 ...
##  $ apg1    : int  8 7 1 8 5 8 9 NA NA 4 ...
##  $ vent    : int  0 0 1 0 1 1 0 NA NA 0 ...
##  $ pneumo  : int  0 0 0 0 1 0 0 NA NA 0 ...
##  $ pda     : int  0 0 0 0 0 0 0 NA NA 0 ...
##  $ cld     : int  0 0 NA 0 0 0 0 NA NA 0 ...
##  $ pvh     : Factor w/ 3 levels "absent","definite",..: NA NA NA NA 2 1 NA NA 1 NA ...
##  $ ivh     : Factor w/ 3 levels "absent","definite",..: NA NA NA NA 2 1 NA NA 1 NA ...
##  $ ipe     : Factor w/ 3 levels "absent","definite",..: NA NA NA NA NA 1 NA NA 1 NA ...
##  $ year    : num  81.5 81.5 81.6 81.6 81.6 ...
##  $ sex     : Factor w/ 2 levels "female","male": 1 1 1 2 1 1 1 NA NA 1 ...
##  $ dead    : int  0 0 1 0 1 0 0 1 1 0 ...

2

Постройте графики плотности распределения для числовых переменных. Удалите выбросы, если таковые имеются. Преобразуйте категориальные переменные в факторы. Для любых двух числовых переменных раскрасьте график по переменной ‘inout’.

df_cleaned2 <- df_cleaned1 %>%
  mutate(
    twn = as.factor(twn),
    vent = as.factor(vent),
    pneumo = as.factor(pneumo),
    pda = as.factor(pda),
    cld = as.factor(cld),
    dead = as.factor(dead),
    apg1 = as.factor(apg1)
  ) %>%
  
  mutate(across(where(is.numeric), ~ ifelse(
    . < quantile(., 0.25, na.rm = TRUE) - 1.5 * IQR(., na.rm = TRUE) | 
    . > quantile(., 0.75, na.rm = TRUE) + 1.5 * IQR(., na.rm = TRUE), 
    NA, .))) %>%
  filter(if_all(everything(), ~ !is.na(.))) %>%
  
  mutate(ID = row_number()) %>%
  relocate(ID, .before = 1) %>%
  
  filter(hospstay >= 0)
  

df2 <- df_cleaned2 %>%
select(where(is.numeric), -ID) %>%
  pivot_longer(cols = everything(), names_to = "variable", values_to = "value")

  ggplot(df2, aes(x = value)) +
  geom_density(fill = "lightblue", alpha = 0.5) +
  facet_wrap(~ variable, scales = "free") +
  labs(title = "Density plots for numeric data", x = "", y = "Density")

plot1 <- ggplot(df_cleaned2, aes(x = pltct, fill = factor(inout))) +
  geom_density(alpha = 0.5) +
  labs(title = "Distribution density for lowph by variable 'inout'", 
       x = "lowph", y = "density", fill = "inout") 
  

plot2 <- ggplot(df_cleaned2, aes(x = bwt, fill = factor(inout))) +
  geom_density(alpha = 0.5) +
  labs(title = "Distribution density for btw by variable 'inout'", 
       x = "bwt", y = "density", fill = "inout") 
  
combined_plot <- ggarrange(plot1, plot2, 
                           ncol = 2,        # Количество столбцов
                           nrow = 1,        # Количество строк
                           common.legend = TRUE,  # Общая легенда для графиков
                           legend = "bottom"      # Размещение легенды внизу
)
combined_plot

df_long <- df_cleaned2 %>%
  select(ID, where(is.numeric)) %>%
  pivot_longer(cols = -ID, names_to = "variable", values_to = "value")

  ggplot(df_long, aes(x = ID, y = value)) +
  geom_point(fill = "lightblue", alpha = 0.5) +
  facet_wrap(~ variable, scales = "free")

3

Проведите тест на сравнение значений колонки ‘lowph’ между группами в переменной inout. Вид статистического теста определите самостоятельно. Визуализируйте результат через библиотеку ‘rstatix’. Как бы вы интерпретировали результат, если бы знали, что более низкое значение lowph ассоциировано с более низкой выживаемостью?

wilcox_test_result <- df_cleaned2 %>%
  wilcox_test(lowph ~ inout) %>%
  mutate(y.position = 7.6) %>%
  print()
## # A tibble: 1 × 8
##   .y.   group1       group2         n1    n2 statistic           p y.position
##   <chr> <chr>        <chr>       <int> <int>     <dbl>       <dbl>      <dbl>
## 1 lowph born at Duke transported   412    71     19947 0.000000957        7.6
ggboxplot(df_cleaned2, x = "inout", y = "lowph", 
          color = "inout", palette = "jco") +
  stat_pvalue_manual(wilcox_test_result, label = "p", tip.length = 0.01) +
  labs(title = "Wilcoxon Rank Sum Test", y = "lowph", x = "inout")

P-value из Wilcoxon теста (5.77e-08) значительно меньше 0.05, следовательно есть статистически значимая разница в значениях lowph между двумя группами: born at Duke и transported. Если более низкие значения lowph ассоциированы с худшей выживаемостью, можно предположить, что младенцы из born at Duke имеют лучшие шансы на выживаемость по сравнению с группой transported.

4

Сделайте новый датафрейм, в котором оставьте только континуальные или ранговые данные, кроме ‘birth’, ‘year’ и ‘exit’. Сделайте корреляционный анализ этих данных. Постройте два любых типа графиков для визуализации корреляций. Постройте иерархическую кластеризацию на этом датафрейме.

df4 <- df_cleaned2 %>%
  select(
    where(is.numeric),
    -c(ID, birth, year, exit)
  )

cor_matrix <- cor(df4, use = "complete.obs")
print(cor_matrix)
##            hospstay      lowph       pltct        bwt        gest
## hospstay  1.0000000 -0.2628111 -0.13751643 -0.5222201 -0.43391064
## lowph    -0.2628111  1.0000000  0.29279438  0.3103002  0.35928694
## pltct    -0.1375164  0.2927944  1.00000000  0.2180920  0.02837645
## bwt      -0.5222201  0.3103002  0.21809200  1.0000000  0.67667411
## gest     -0.4339106  0.3592869  0.02837645  0.6766741  1.00000000
plot4.1 <- df4 %>% 
  cor() %>% 
  corrplot(
    order = 'hclust'
  )

plot4.2 <- df4 %>% 
  cor() %>% 
  network_plot(min_cor = .0)

print(plot4.1)
## $corr
##            hospstay       pltct      lowph        bwt        gest
## hospstay  1.0000000 -0.13751643 -0.2628111 -0.5222201 -0.43391064
## pltct    -0.1375164  1.00000000  0.2927944  0.2180920  0.02837645
## lowph    -0.2628111  0.29279438  1.0000000  0.3103002  0.35928694
## bwt      -0.5222201  0.21809200  0.3103002  1.0000000  0.67667411
## gest     -0.4339106  0.02837645  0.3592869  0.6766741  1.00000000
## 
## $corrPos
##       xName    yName x y        corr
## 1  hospstay hospstay 1 5  1.00000000
## 2  hospstay    pltct 1 4 -0.13751643
## 3  hospstay    lowph 1 3 -0.26281113
## 4  hospstay      bwt 1 2 -0.52222009
## 5  hospstay     gest 1 1 -0.43391064
## 6     pltct hospstay 2 5 -0.13751643
## 7     pltct    pltct 2 4  1.00000000
## 8     pltct    lowph 2 3  0.29279438
## 9     pltct      bwt 2 2  0.21809200
## 10    pltct     gest 2 1  0.02837645
## 11    lowph hospstay 3 5 -0.26281113
## 12    lowph    pltct 3 4  0.29279438
## 13    lowph    lowph 3 3  1.00000000
## 14    lowph      bwt 3 2  0.31030019
## 15    lowph     gest 3 1  0.35928694
## 16      bwt hospstay 4 5 -0.52222009
## 17      bwt    pltct 4 4  0.21809200
## 18      bwt    lowph 4 3  0.31030019
## 19      bwt      bwt 4 2  1.00000000
## 20      bwt     gest 4 1  0.67667411
## 21     gest hospstay 5 5 -0.43391064
## 22     gest    pltct 5 4  0.02837645
## 23     gest    lowph 5 3  0.35928694
## 24     gest      bwt 5 2  0.67667411
## 25     gest     gest 5 1  1.00000000
## 
## $arg
## $arg$type
## [1] "full"
print(plot4.2)

5

Постройте иерархическую кластеризацию на этом датафрейме.

df_scaled <- scale(df4)
  rownames(df_scaled) <- df_cleaned2$ID
dist_matrix <- dist(df_scaled, method = "euclidean")
hc <- hclust(dist_matrix, method = "ward.D2")

fviz_dend(hc, 
          cex = 0.6)
## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## ℹ The deprecated feature was likely used in the factoextra package.
##   Please report the issue at <https://github.com/kassambara/factoextra/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

#fviz_dend(hc, 
          #k = 3, 
          #k_colors = "jco",
          #type = "phylogenic", 
          #repel = TRUE) # Избежать наслоения лейблов

6

Сделайте одновременный график heatmap и иерархической кластеризации. Интерпретируйте результат.

library(pheatmap)
pheatmap(df_scaled,
         clustering_method = "ward.D2",  # Метод кластеризации
         scale = "none",                # Данные уже нормализованы
         color = colorRampPalette(c("navy", "white", "firebrick"))(50),  # Цветовая палитра
         main = "Heatmap with Hierarchical Clustering")

Субъекты с преимущественно низкими значениями hospstay могут представлять наблюдения с большими значениями bwt, gest и lowph и наоборот.

7

Проведите PCA анализ на этих данных. Проинтерпретируйте результат. Нужно ли применять шкалирование для этих данных перед проведением PCA?

df.pca <- prcomp(df4, 
                scale = T)
summary(df.pca)
## Importance of components:
##                           PC1    PC2    PC3    PC4    PC5
## Standard deviation     1.5461 1.0294 0.8393 0.7541 0.5263
## Proportion of Variance 0.4781 0.2119 0.1409 0.1137 0.0554
## Cumulative Proportion  0.4781 0.6900 0.8309 0.9446 1.0000
fviz_eig(df.pca, 
         addlabels = T, 
         ylim = c(0, 50))

fviz_pca_var(df.pca, col.var = "contrib")

Первая компонента объясняет 47.8% дисперсии данных, первые две компоненты объясняют 69% дисперсии данных, можно сделать вывод что их можно использовать для дальнейшего анализа. Шкалирование применять нужно, так как данные имеют разные единицы измерения и разную дисперсию.

8

Постройте biplot график для PCA. Раскрасьте его по значению колонки ‘dead’.

library(ggbiplot)

biplot <- ggbiplot(df.pca, 
         scale=0, alpha = 0.7,
         groups = as.factor(df_cleaned2$dead))+
  scale_color_discrete(
    name = "Status",
    labels = c("0" = "Alive", "1" = "Dead")
  ) +
  labs(title = "PCA Biplot colored by 'dead'")
biplot

9

Переведите последний график в ‘plotly’. При наведении на точку нужно, чтобы отображалось id пациента.

library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
pca_data <- as.data.frame(df.pca$x)  # Извлекаем главные компоненты
pca_data$ID <- rownames(df_cleaned2)         # Добавляем ID пациентов
pca_data$dead <- as.factor(df_cleaned2$dead)
pca_data$dead <- factor(pca_data$dead, levels = c("0", "1"), labels = c("Alive", "Dead"))

plot9 <- ggplot(pca_data, aes(x = PC1, y = PC2, color = dead, text = ID)) +
  geom_point(size = 2, alpha = 0.7) +
  labs(title = "PCA Biplot colored by 'dead'", 
       x = "PC1", 
       y = "PC2",
       )+
scale_color_discrete(
    name = "Status",
    labels = c("0" = "Alive", "1" = "Dead"))

# Преобразование в интерактивный график с помощью plotly
plotly_plot <- ggplotly(plot9, tooltip = "text")  # 'text' отображает ID пациента
plotly_plot

10

Дайте содержательную интерпретацию PCA анализу. Почему использовать колонку ‘dead’ для выводов об ассоциации с выживаемостью некорректно?

Вообще (PCA) используется для уменьшения размерности данных, выявления главных направлений вариации в данных и для упрощения визуализации данных. Как было сказано ранее, первая компонента объясняет 47.8% дисперсии данных, первые две компоненты объясняют 69% дисперсии данных, что является хорошим показателем для дальнейшего анализа. Самые длинные стрелки у переменных hospstay и gest, они вкладывают больше всего информации в объяснении вариации, также видно что они направлены в одну сторону, что свидетельствует об их корреляции. Вклад hosptstay, bwt и gest сильнее в компоненту PC1, pltct и lowph влияют на PC2. Из графика видно, что точки с Alive и Dead значениями не образуют кластеров, что указывает отсутствие тенденции разделения выживших и умерших на основе первых двух главных компонент.

PCA основывается исключительно на входных переменных (hosptstay, bwt, gest, lowph, pltct) и их вариации. Переменная dead используется только для визуализации и никак не влияет на формирование главных компонент. Чтобы проверить, связаны ли переменные с выживаемостью корректнее использовать модели, такие как логистическая регрессия.